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We study the statistical mechanics of a class of problems whose phase space is the set of per- 
mutations of an ensemble of quenched random positions. Specific examples analyzed are the finite 
temperature traveling salesman problem on several different domains and various problems in one 
dimension such as the so called descent problem. We first motivate our method by analyzing these 
problems using the annealed approximation, then the limit of a large number of points we develop 
a formalism to carry out the quenched calculation. This formalism does not require the replica 
method and its predictions are found to agree with Monte Carlo simulations. In addition our 
method reproduces an exact mathematical result for the Maximum traveling salesman problem in 
two dimensions and suggests its generalization to higher dimensions. The general approach may 
provide an alternative method to study certain systems with quenched disorder. 



I. INTRODUCTION 



The statistical mechanical approach to the study of optimization problems has lead to progress in a number of 
ways. The approach is based on identifying the cost function, which needs to be minimized, with the energy of a 
physical system whose phase space is equivalent to the free adjustable parameters in the optimization problem. The 
zero temperature energy of the resulting physical system thus corresponds to the optimal solution. This formulation 
can be exploited in two ways. First, physically motivated minimization techniques such as simulated annealing can 
be applied to optimization problems 1], often leading to near optimal solutions. Secondly the statistical mechanical 
approach can also be used to carry out computations of average or typical values of optimal solutions, where the 
non-adjustable parameters (describing the realization of the instance) in the system are taken to be quenched random 
variables Q- The replica and cavity methods, which are much used in the theory of spin glasses, have been successfully 
exploited to study statistical properties in wide range of optimization problems [3, 0, 0,0 Often optimization problems 
have a phase space which is equivalent to permutations or partitions of the integers and these problems are referred 
to as combinatorial optimization problems. One of the most famous of these combinatorial problems is the traveling 
salesman problem (TSP). Here the problem is to find the minimal circuit length to visit N cities or points where 
the distance between the points i and j is given by dij. The order in which the cites are visited is encoded in a 
permutation a 6 S^r where £jv is the group of permutations of N objects. For a given permutation 

D{a) = Y,d tTi , rTi+1 , (1) 

i 

is the corresponding total distance traveled. When the dijS are chosen from some quenched distribution the problem 
is referred to as the stochastic TSP. The most natural form of the TSP is the Euclidean TSP where the cities are 
points ri,r2---rjv in some connected domain T> in R d and each point is independently distributed from the others 
with the same probability density function p q (r). The distance between the points i and j is simply the Euclidean 
distance on M. d given by dij = |r, — Tj\. It was shown 6] that for N — > 00 the minimal path Dm behaves as 

THTr-jS(d) / d d r (p^r)) 1 "* (2) 

with probability one. Here (3(d) is a constant depending on the dimension of the space d but independent of p q . The 
stochastic TSP has also been studied under the random link hypothesis where the dij are all uncorrelated (up to any 
symmetry requirement). Clearly in this random link version the triangle inequality is not respected. This version 
has been intensively studied 0, H H E3, EH 

and its analysis is greatly simplified by the lack of correlation between 
the dij which makes the taking of the disorder average quite straight forward. A, somewhat perverse, variant of the 
TSP is one where one asks for the maximal tour, this is called the Maximum TSP ,12] and is, for obvious reasons, 
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FIG. 1: Ring polymer on a two dimensional substrate where each monomer (filled circles) is attached to an impurity (crosses) 



sometimes referred to as the taxicab rip off. In the statistical mechanical formulation if one looks for the maximal 
tour, one keeps the same cost function but changes the sign of the temperature. 

In the class of problems we shall study in this paper, N points {ri,r2 • • • , r^v} are chosen independently in some 
domain T> S M. d with probability density p q (j)- Again the dynamical phase space for the problem is taken to be all 
permutations of the order of these points, a € Sat. The Hamiltonian for the system is defined to be 



N 



%) = ) /r„,- Vl , (3) 



Cyclic boundary conditions ro = are imposed. In this context the Euclidean TSP problem corresponds to the zero 
temperature limit of the case where V(r) = r. 

A physical realization of the system is one where the are impurities where the monomers of a polymer loop are 
pinned, and only one monomer can be pinned per impurity. The potential V represents effective interaction between 
neighboring monomers on the chain. For instance V(r) = Ar 2 /2 corresponds to the Rouse model of a polymer 
chain[13j. For instance in Fig. we represent the system of a polymer on a two-dimensional substrate where the 
monomers, shown as filled circles, attach themselves to pinning sites (shown as crosses). 

The canonical partition function for these problems is given by the following sum over all permutations 

Zn = W\ £ ex P(-/3^))- ( 4 ) 

' crGSjv 

Since the number of permutations grows as N\ the entropy is non-extensive and behaves as JVlniV, but here we insert 
a factor of 1/iV! to absorb it. 

To compute the average energy per site it is necessary to work out the quenched free energy of the system 

1: 



F N = --hy(Z N ), (5) 

where the over-line denotes averaging with respect to the quenched joint probability density function of the random 
sites rj. The energy per site is then evaluated as 

Our method will be shown to be exact in the limit of large N while keeping the domain T> fixed. However this 
large N scaling is not the one needed to obtain the quantity (3{d) in Eq. (J2J. When the probability density p q (r) is 
flat and V(r) is an attractive potential, as is the case for the ordinary TSP, the minimal energy configuration is one 
where links are always of the order of the minimal separation between points, that is to say 0(p^ 1 ^ d ) where p is the 
density of points, and thus the ground energy per site is of the order 

ecs » Vip- 1 /*). (7) 
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If one uses the TSP potential, V^r) = |r|, in the above, one recovers the scaling of Eq. (J2J. In these attractive cases 
the ground state energy per site is zero in the thermodynamic limit and in order to extract an extensive result the 
energy of the system must be scaled appropriately with TV 3J. If the potential V is repulsive, as in the Maximum 
TSP, then links will be typically of the domain size T> and if this domain size is O(l) (in the sense that it does not 
scale with N) then the ground state energy per site will be O(l) without the need for any special scaling. 

We emphasize that our method is exact in any dimension, however many of the examples we give will be in one 
dimension, where many explicit results can be obtained. In one dimension we remark that any attractive choice of the 
potential V corresponds to a choice of cost function for the computational problem of sorting random data elements 
into increasing order. The performance of local physical Monte Carlo algorithms has been analyzed in problems with 
these cost functions |16| . demonstrating the pitfalls in using such algorithms to search for the optimum. Indeed tree 
based sorting algorithms are much more efficient |15| . 

The main potential we shall investigate is 

V(r) = |r|, (8) 

and, as mentioned previously, this potential is of particular interest as it arises naturally in the Euclidean TSP as its 
ground state at positive temperature is the shortest circuit visiting each of the points once and only once. We will 
discuss this problem in detail in both one and two dimensions on several domains with different topologies. Some of 
these domains are not simply connected and besides the average energy observable, that is our main focus, we show 
how to compute the statistics of the winding number of the path around the domain. Of course the solution to the one 
dimensional TSP is obvious, one starts with the leftmost point and works along to the rightmost, giving an average 
ground state energy per site of ecs — 0. Using the scaling of this paper the ground state energy of the ordinary TSP 
per site is cqs = in all dimensions, as can be seen from Eq.(J2J). At negative temperature (or where the sign of V 
is inverted) the corresponding ground state corresponds to the solution of the Maximum TSP where one requires the 
maximal distance taken to complete a circuit visiting all the points once and only once. The solution here is not quite 
so obvious but we shall see that ecs = 1/2 is the average ground state energy per site. For a uniform distribution of 
points the interested reader may verify that this average value of ecs may be achieved by a greedy algorithm, which 
starts at the leftmost point, goes to the rightmost point, returns to the next leftmost point and so on. We emphasize 
however that the main point of this paper is to solve the finite temperature statistical mechanics of these models at 
all temperatures. 

We shall also consider the descent problem 0] as a fully soluble system in which all our equations can be resolved 
analytically. This is a one dimensional system described by the potential 

V{x) = 6{-x). (9) 

for x G [0,1]. In fact the original formulation of the descent model, is equivalent to one where x% are chosen to 
deterministically as Xi = i/N. However because of the scale free nature of the potential V, all models with an 
arbitrary continuous distribution of the Xi are in fact equivalent. The energy of the permutation is thus the number 
of points where x ai is greater than the point x ai+1 , the point which follows it on the polymer ring. The ground state 
of the system is simply the permutation in which the x ai appear in increasing order and has corresponding energy 
per site ecs — 0. 

To further test the validity of our method, we have considered harmonic potentials V(r) = |r| 2 in various dimensions. 
In one dimension with domain Xi € [0, 1], we have also studied the potential V(x) = — ln(|x|) at positive temperatures 
(where the model is defined). 

In all cases studied the analytic predictions were confirmed by Monte Carlo simulations and in some cases by 
extrapolating the results of exact enumeration for systems of small size. In addition the behavior of ecs for repulsive 
potentials V in one dimension is analyzed via a zero temperature analysis, the results are discussed in terms of the 
corresponding optimal paths. We also examine the Maximum TSP in higher dimensions and show that we recover an 
exact mathematical result for the average length of the optimal path in two dimensions. We are able to use our method 
to predict the corresponding optimal path length in higher dimensions. A fascinating aspect of this analysis is that 
in addition to providing the average optimal path length, the saddle point equations we derive in the thermodynamic 
limit seem to suggest the heuristic one should use to search for the optimal path. 

The rest of paper is organized as follows: in section (II) we study the annealed approximation in order to motivate 
the new method for the quenched case which is presented in section (III). Section (IV) is concerned with the zero 
temperature limit and results for the average length (or ground state energy) of optimal paths. 

A brief description of our method has appeared in 01 an d a comment on the technique can be found on 
:/ /jc-cond-mat. bel-labs.com 
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II. THE ANNEALED APPROXIMATION 
A. General formalism 

Here we shall analyze the statistical mechanics of this general class of problems in the annealed approximation 
which amounts to setting 

F N ^FST n = ~]n(Z N ), (10) 

This approximation is in general doomed to failure for the following reason. In the annealed approximation the 
quenched variables are no longer quenched and will evolve dynamically in order to decrease the free energy of the 
system. The configuration of variables which dominates the thermodynamics will generically be atypical of the initial 
quenched distribution and will usually be of measure zero. For instance, if we consider the one-dimensional TSP at 
negative temperature it is clear that the maximal circuit, averaging over all permutations and positions Xi, will be one 
where half the Xi are at the point x — and the other half at x = 1 . This will allow a ground state energy of ecs = 1 
from applying the greedy algorithm mentioned in the Introduction. However this configuration is of measure zero if 
the quenched distribution is uniform on [0,1]. Despite this deficiency, the formalism below will have an important 
bearing on our subsequent development of the quenched calculation and indeed its physical interpretation. 

The partition function is averaged over all rj, i = 1 . . . N, with periodic boundary conditions, i.e. ro = rjv- Upon 
this averaging all permutations become equivalent and we obtain 

N / N \ 



Z N = / Y[d d n cxp -jS^Vfrrrj-i) • (11) 



i=i 



The above averaged partition function can be evaluated using standard transfer operator techniques: 

= TrT N = J d d rT N {r,r), (12) 

where T is the operator 

T(r,r')=cxp(-/?V(r-r')). (13) 
In the limit of large N, taking the size of the domain to be normalized to unity, we find 

Z^=^a, (14) 

where A a is the largest eigenvalue of the operator T (we use the subscript a to indicate a quantity evaluated in the 
annealed approximation) . The corresponding right and left eigenfunctions \ of T obey 

/i a) (r) = A" 1 J d d r> cxp (-f3V(r r')) if V) (15) 

/| a) (r) = A" 1 J d d r> exp {-(3V{r' rj) ft\r>) (16) 

These eigenfunctions are identical for a symmetric potential and we chose the eigenfunctions to be normalized. 

In this cigen-system, A a can have several solutions but that with the maximum value of A a dominates the partition 
function at large N . Moreover the eigen- functions f^ L corresponding to the largest eigenvalue must be positive by 
the Perron-Frobenius theorem. The annealed approximation for the energy per site is consequently obtained as 

_ 01n(Aq) (m 

The annealed density of points on the path at the point r is given by 

1 N 
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where the angled brackets indicates the Gibbs ensemble average over the annealed points. From the periodic boundary 
conditions, all points are equivalent and we have 

T N (r r) 

Va {v) = <<5(r-n)) - 



Zn 

= f ( R a) (r)f { L a \r), (19) 

where again we have taken the thermodynamic limit. 

In this annealed approximation the probability density of the points is thus given by Eq. I|19f) . I n general we will 
find that 

Pa{v)^ Pq {v) (20) 

as the variables evolve dynamically. 

A general expression for the energy at high temperature can be obtained by expanding the averaged partition 
function for small j3. The first two terms of this expansion are: 

e a = J d d rd d r' V(r-r')+{3 3 (^j d d rd d r' V(r - r')\ -2 J d d rd d r'd d r" V(r - r')V(r' - r") - J d d r'd d r V 2 (r - r') 
A slightly more involved calculation for the quenched case yields the differing expansion: 

e = J d d rd d r' V(r-r')+/3 2 J d d rd d r'd d r" V(r - r')V(r' - r") - J d d rd d r' V 2 (r - r') - (J d d rd d r' V(r - r') \ 

( 22_ ) 

Thus, in general, there will be a difference between the annealed approximation and quenched result at any finite 
temperature. We will later check our method for the quenched case by seeing that it reproduces the second form Eq. 

B. The descent model 

We start by considering the annealed approximation for the descent model. Differentiating the first of Eq. I|16() 
with respect to x yields 

elf a 

^ = A- 1 /^(l-exp(- y 5)) (23) 

This has the solution 

/£ = Cexp(Kz) (24) 

where k = A a ; 1 (l — exp(— /?)) and C is a constant of normalization. This solution is then substituted into the original 
integral equation to yield A a = (1 — exp(— /?))//?. Then Eq. I|17|l gives the annealed energy to be 

(25) 



(3 exp(/3) - 1 

The result Eq. (|25|l is in fact identical to that obtained from the exact solution to the descent problem, obtained 
via more lengthy combinatorial methods 14] . The annealed approximation thus leads to the exact energy per site 
for this problem. This exactness is straightforward to understand. In this problem there is no length scale in the 
potential V and only the order of the points determines the energy, clearly one would obtain the same energy from 
any continuous distribution of points selected independently with density p q (x) on [0, 1]. We note that the solution 
for the other (left) eigenvector /£ is 

ft = C cxp(-Kx) (26) 

and we obtain find p a [x) = 1 = p q (x), thus the annealed distribution agrees with the quenched one. This is an 
autoconsistency of the annealed approximation which leads to it being exact. 
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C. The One-Dimensional TSP 



In this section we restrict ourselves to one dimension and consider the unit interval with x £ [0, 1]. We note that 
the potential V is symmetric and thus may write f^ 1 ' = = f^ a \ 
We proceed by differentiating Eq. twice to obtain 

^_ /3 2 / (a) +2 £ / (a) =0i (2?) 

which has solution 

/(») = C (exp(ujx) + Acxp(-cjx)) (28) 



with 



Now, A a must be positive, so we have 



lj 2 < [3 2 for (3 > (30) 
lj 2 > (3 2 for /3 < (31) 

Substituting the solution Eq. I|28|) back into Eq. I|16|) we find the condition 

ex P (2.) = (^) 2 . (32) 

We note that a solution of Eq. I|32l) is a; = 0. However the corresponding solution for f( a > would be of the form 
f^(x) = Ax + B. From above one must also have that df^ /dx is continuous and the clear symmetry f^{x) = 
/( a )(l — x) means that A = in this solution. One can verify that f^(x) = B is not a solution. Hence the solution 
must have uj 0. For (3 > one finds that to — i(3z where z is the smallest positive solution of 

z = cot(/3z/2), (33) 

and shows no discontinuities as (3 varies. The annealed energy per site is then given by 

1 z 2 

(34) 



a S3 (l + i/3(l + 2 2 )) 
For (3 < we find that ui = (3z where z is the root of 

z = coth(-/3z/2). (35) 



and here the annealed energy per site is 



In both cases, as (3 — > we find that 



1 z 2 

J3 + (1+1/3(1-^)) (36) 



6a(0) - i (37) 

which agrees with the exact high temperature result Eq. (|21|l . 

In the positive temperature case the smallest positive solution z* to Eq. I)33|l is such that z* £ (0, 7r//3) and thus 
z* — > as [3 — > oo and thus e a rs — ► 0; this clearly agrees with the limiting behavior of the corresponding 
quenched case. In the case of negative temperature we find that the solution of Eq. (|33|l behaves as z* — > 1 (plus an 
exponentially decaying correction), thus yielding z* sa 1 + 1//3 — * 1 which implies e a — » 1. This latter result is clearly 
not the correct result for the quenched case as the distribution of the Xi has evolved to permit this maximal energy 
configuration as explained in the Introduction. Indeed we find from Eq. (|28|l that p a (x) = f 2 ( x ) 7^ 1 and becomes 
peaked at the boundaries. 
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D. The TSP on a Ring 



When the domain of the one-dimensional TSP is a periodic ring rather than the unit interval with boundaries 
considered in the last section, the analysis becomes simpler, moreover there is an interesting new observable, the 
winding number. For a ring domain, the shortest path is sometimes the other way round the ring, so the potential is 

^> = {'H4 EE>I (38> 

Taking proper account of the contributions from discontinuities in the derivative of this potential, the equivalent 
formula to (|27|l gains a term on the right hand side: 

_ 0> f (a) {x) + 2 £/W(x) = (a; + 1/2), (39) 

where the eigenfunction is now a periodic function with period 1. In this case, in contrast to the situation with 
boundaries, the original integral equation (|16J) admits a constant solution f a = l. For positive {3, this must correspond 
to the largest eigenvalue which is: 

2(1 - e- / 2 ) 

K . ?<1_£_> ,40, 

The annealed energy per site is then given by 

£a = p ~ 2(e/3/ 2 - 1) (41) 

The high temperature limit of e a now takes the value 1/4. 

Notice that the constant solution indicates that the annealed density of points is the same as that of the desired 
quenched distribution, so the annealed approximation is exact in this case. This result follows from the symmetry 
of the domain and continues to hold for certain other closed domains in higher dimensions. The higher dimensional 
cases will be treated in the later section on quenched models. 

A new observable, the winding number, arises for this domain. The winding number counts the number of times a 
particular path goes around the ring and can be written as the sum of contributions from each step between points 
in the same way as the original Hamiltonian. 

JV 

W(a)=J2w(x cri -x lTi _ 1 ), (42) 
i=l 

Where: 

{l + x, for x < — |; 
x, for \x\ < |; (43) 

— l + x, for x > |. 

In fact, since the cyclic boundaries lead to — x (Ti _ 1 ) = 0, the x part of the contributions can be dropped and 

the winding number may be written: 

N 

W{a) = Y,-° ( x °* - x "*-i - V2) + (-z<r 4 + x ai _, - 1/2) , (44) 

8=1 

In this form it shows some similarity with the potential for the descent model. In particular, the asymmetry will make 
the right and left eigenfunctions differ. 

We compute the expectation values of the winding number by taking the Boltzmann weight of the configuration a 
to be 



exp(-(3H(cr) -jW(a)) 



(45) 
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FIG. 2: Expectation value for the fluctuations in winding number per site, {W 2 )/N, for the one dimensional TSP on the ring 
domain as a function of f3 (dotted line) compared with the Monte Carlo simulations (crosses). Negative ft corresponds to the 
Maximum TSP problem. 

The analysis developed above holds for this slightly more general case and we will obtain expectations of the winding 
number by differentiating the partition function with respect to 7 before setting 7 to zero. 

By considering the equations obtained by differentiating the eigenvalue equation it becomes clear that a solution of 
the form: 



f% = Cexp(ra) 



(46) 



should be sought. Indeed, this form is a solution of the integral equation provided k = 7. As was the case for the 
descent model, the left eigenvalue is of the form /£ = C"exp(— nx), so the annealed density becomes constant and 
this solution is thus also valid for the quenched problem. The corresponding eigenvalue is given by 



A = 



2ft 



,-/3/2 



ft 2 - 7 2 ft 2 



r 



((/3 + 7 )^ /2 + (/?-7)e^ /2 ) 



(47) 



which correctly reduces to Eq. 1)41) [) when 7 is set to zero. 

As should be exppected, the mean value of the winding number vanishes since there is nothing that prefers winding in 
one direction over the other. Even in the low temperature limit, there is no symmetry breaking in this one-dimensional 
system. On the other hand, the fluctuations provide a non vanishing observable: 



__ = N_ (8e^ -8- Aft -ft 2 ) 



Aft 2 



( e /3/2 



1 



(48) 



In the high temperature limit the fluctuations per site become 1/12 and for large negative ft the limit is 1/4. 

We have performed Monte Carlo simulations to test this prediction using the quenched model which we have already 
argued has the same value of observables. These results are show good agreement as displayed in Fig. J2J. 



III. THE QUENCHED CALCULATION 
A. General formalism 



The order of the points is unimportant for the statistical mechanics of this problem because the phase space is 
all their possible orderings. The relevant disorder is thus clearly fully determined by the A*" position vectors of the 
sites. These positions are encoded in the, unaveraged, density of points in space 



1 N 



(49) 
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where we have used the superscript q above to emphasize that the points rf^ are quenched. We note that by definition 
we have p q {r) = ~p q {v). and in the limit of large A we expect that 

' d d r p q (r)h(r) = I d d r p q {v)h{r) + 0(1/ VF) (50) 



for suitably well behaved functions h. If is the site visited by the polymer at step i in a system which has G T> 
then the partition function of the permutation problem can be written as 



' ^ i=l r r \ i ) V i / 



(51) 



The above can be derived by considering a discrete version of the problem where one has n(j) sites at the points 
r(j). A path is specified by the possible sequences ri, T2, • • • rjy; however at the visit to the site r(j) there are n(j) 
possible points to chose from and thus any path has a degeneracy Y[j n(J)\ in order to have the same phase space 
as the permutation problem. In addition each site r(j) can only be visited n(j) times, explaining the delta function 
constraint above. Another way of obtaining the factor in Eq. <|51|) is to note that Zm{0) = Cn£n((3) where 



En = /n d d r i n*^ 9 (r)-E*( r - r *)) exp^-Z^F^-r^ 



(52) 



is the partition function of the system up to a temperature independent entropy/degeneracy contribution Cat. Clearly, 
as defined here, Z(0) = 1 which implies Z((3) = S(/3)/H(0). We shall later that, in the limit of large A we have 



Sjv(O) = exp f—N J d d r p q {v)[\n{p q (r)) - I 



(53) 



which in the large N limit, via Stirling's formula, recovers Eq. (|51|l . 

The partition function S^r may be written using a Fourier representation of the functional constraint: 



E N = J d[n\ exp (n J d d i 



r/i(r)p fl (r)j Z N , (54) 

where each integration over /i(x) is up the imaginary axis. The object Zn is similar to the annealed partition function 
considered in the previous section, but with an r dependent chemical potential. It is defined as: 



„ N / N N \ 

2n = / n ddr * ex p -p v ^ - r '-i) - E ^ ■ ( 55 ) 

J i=l \ i=\ i=0 I 



For large N we may use the relation Eq. H5l)fl . neglect the terms 0(vN), and then the partition function in Eq. I|54l) 
can be evaluated by the saddle point method in the limit where N — > oo keeping T> fixed. The saddle point equation 
is 

1 5\uZ N 1 , .. , . <51n(A 9 ) 

^ (r) = "aT^T = A ( g^- ri)> = Mr) - (56) 

where the above expectation, is in the system with partition function Z^ defined in Eq. (|55|l . 

Physically this approach can be thought of as choosing a site dependent chemical potential fi which fixes the density 
of the annealed calculation to be the same as that of the quenched one p q (r) — £> (r). This idea was used sometime 
ago in an approximative sense where low order moments, not the whole distribution, were fixed in this way [TsUIi^. 

To proceed, we find an expression for p a (r) using similar techniques to those employed in the last section. We find 



l N = TrT N (57) 



where 



T(r, r') = exp (-/x(r)/2) exp {-f3V{r - r')) exp (-/*(r')/2) 



(58) 
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The right and left ground state eigenfunction, corresponding to the maximal eigenvalue A^, obey 

/£>(r) = \- l exp(-»(T)/2) J d d r' exp (-0V(t - r')) exp (-/x(r')/2) /^(r') (59) 

/| 9) (r) = X-'expi-fxir)^) J d d r' exp (-/3y(r' - r)) exp (-//(r')/2) /| ?) (r') (60) 
and by a similar calculation to that of the annealed case we have: 

p,(r) = Po (r) = -^[^ = /W(r)/«(r). (61) 

In the case of a symmetric potential where V(r) = V(— r), we have that — = and thus = y/p^, note 
this ensures that is the eigenfunction corresponding to the maximal eigenvalue as we note that it is positive and 
then appeal to the Perron-Frobenius theorem. Substituting the above into the saddle point equation gives 

^) = A- 1 cxp(-^) j d d r' expH3V(r-r'))cxp(-^) v /^0 (62) 

we can thus write exp(— = ^/p q {v) / s\ q (r) where s\ q (r) obeys 

sx g (r) = A" 1 / d d r' exp (-0V(t - r')) ^\ (63) 

Substituting this back into the action we obtain 
HZn) 



N 



2 / d d r Pq (v) In (s Xq (r)) + ln(A ? ) - / d d r Pq (r) In(p 9 (r)), (64) 



the last term which is independent of explains the presence of the combinatorial term in Eq. I|51|l . However from 
Eq. H63fl w e see that there is a whole family of solutions {sA 5 (r), X q } which are related by s\ q = a 1 / 2 s a \ ql for a > 
and in addition these solutions all have the same action. This apparent zero mode is an artifact introduced by the 
fact that the constraint N j d d r p q = J d d r 6(r — r.;) is automatically satisfied. Thus we may chose X q = 1. In the 
case of a uniform distribution on a domain of unit volume this leads to our final result 



-2± 
dp 



d d r ln(s(r)) 



= /^dV^-^-^-^ (65) 
J s(r)s(r') 

where s obeys 

J s(r') 

and we have specialized to the usual case when the distribution of quenched points is uniform. Recall that we have 
set the size of the domain to one, and some scaling is needed when this is not the case. 

It is possible to check our result via direct comparison with the high temperature expansion Eq. (|22J) written down 
at the end of section (II). Equation l|66J) can be solved perturbatively as a power series in f3 by writing 

s(r) = l + /3 Sl (r)+/3 2 s 2 (r)--- (67) 

Substituting this expansion in to Eq. (|66|l yields 

s x (r) = - J d d r' V(r - r') + i J d d r'd d r" V(r" - r') 

s 2 (r) = ~ J d d r' V 2 (r-r')+ J d d r' Sl (r')V(r - r') (68) 

- j J d d r'dz V 2 {r' J d d r'd d r" Sl (r')V(r' ~ r") + 1 / d d r' si(r') (69) 

To order (3 Eq. 1651) then yields 

e = -2 J d d r L (r) + 2/?(* a (r) - ^) + 0(/3 2 ) (70) 
and substituting Eqs. I|69|) and l|69|l in the above we recover the result Eq. 122(1. 
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Rouse Polymer in 2-Dimensions 






Beta 



FIG. 3: Theoretical prediction for the average energy e for two dimensional quenched polymer model as a function of /3 (solid 
line) compared with the Monte Carlo simulations (solid circles). Negative /3 corresponds to the maximal problem. Error bars 
based on 20 realizations of the quenched points are smaller than the symbol sizes. Here the domain T> is unbounded but p q is 
Gaussian and centered at the origin 

B. Quenched Rouse Polymer 

Before moving on to consider the TSP, we first analyze a model with a harmonic potential to demonstrate an 
analytic solution of the quenched equations. This system, in the annealed case, is used to model a polymer |13|. To 
avoid boundaries, which prevent an analytic solution, we treat the case where the quenched distribution of points is 



Gaussian: p q {r) 



-r-72 



/(2tt)s and we work in arbitrary dimension d. The quenched equation for s(r) becomes: 



s(r) 



d d r 
(2^1 



,-r' 2 /2 p 



exp (— /3|r — r' 



s(r') 



So we search for a Gaussian solution, 



This is satisfied provided: 



s(r) 



1 = 



I ( i + 2/3 - ^TTW 

-d/4 



2 

2 d / 4 (l + Vl + 4/3 2 ) 



Inserting these into the saddle point action we obtain the average energy per site as: 



= d 1 - 



2/3 



2/3 



y/1 + 4/? 2 v/l + 4/3 2 (l + v/1 + 4/3 2 ) 



(71) 

(72) 

(73) 
(74) 

(75) 



This has the correct d/2/3 behavior at large /3 and the value e = d at infinite temperature, as can be checked directly 
by a Gaussian average. In contrast to the situation for directed polymers, there is no evidence for a phase transition 
in this quenched case. 

In figure @ we show the average energy in the two-dimensional case and compare it with Monte Carlo simulations. 
The agreement is good, though it should be noted that at the edges of the plot, for large |/3|, long runs (10's of millions 
of steps) with large N are required to see accurate agreement. 



TSP in One and Two Dimensions 



In this section we consider predictions for the TSP in both one and two dimensions. In view of the nonlinear nature 
of equations I|b5fl an d l|bb(l we have not found any non-trivial analytic solutions for the TSP potential. Our primary 
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tool is the iterative numerical solution of Eq. I|66|l which is stable and can be solved to any required accuracy. However, 
for closed symmetric domains (we shall consider a one dimensional ring and a disc and torus in two dimensions), a 
constant solution exists and some analytic progress is possible. To see this, we simply require a domain such that 
that the origin of the integration in Eq. Ijtitjf) can be shifted to yield: 



The significance of this observation is that the annealed approximation is exact for these domains. Indeed this equation 
is exactly the annealed eigenvalue equation for a constant eigenfunction and the relationship between the value of the 
constant s and the eigenvalue A of the annealed approximation is simply A = s 2 . In all cases the normalization is 
A(/3 = 0) = 1. This conclusion is consistent with what is known from analysis of the independent link version of the 
TSP. Evidently, the geometry of independent links has no boundaries and the analysis of this problem in the same 
scaling limit we consider here also shows that that the annealed approximation is exact 0. 

The expression (|76|) allows a general large (3 expansion for these closed domains. Provided the potential is convex, 
then in this limit it is apparent that the only contribution to the integral is from nearby points. When the potential 
is smooth we can expand to find A ~ e -^°> (/3V")~ d/2 , so the energy is e ~ ^(0) + d/2(3. This expansion is invalid 
for the TSP (and descent model) since the potential is not smooth, and the correct energy is e ~ d/f3. These results 
indicate that in this limit the topology of the domain becomes unimportant and only the dimension is relevant as 
is indeed observed in all the examples below. A similar limit for the maximal problem is dominated by a term 
corresponding to the largest distance two points can be apart from each other, which is sensitively dependent on the 
topology of the domain. 



For the case of the ring domain, constant s is a solution and as expected this reproduces the energy obtained via 
the annealed calculation (|41|) . 

For the unit interval, we use Eq. I|65|l to predict the energy by iteratively solving Eq. (|66J) . 

To test these predictions we have carried out Monte Carlo simulations of the TSP problem for system sizes of 
N = 5000 and compared the average energy measured after equilibrating the system over 5000 Monte Carlo sweeps 
and measuring the average energy over a subsequent 5000 Monte Carlo sweeps. The disorder average was carried 
out by averaging the results over 20 independent realizations of the disorder. The standard move was taken to be 
a random transposition of a pair of points in the permutation and the acceptance of the move was chosen with the 
Metropolis rule. The results for both the ring and line domains are shown in Fig. Q compared with the predictions. 
We see that for all temperatures the agreement is excellent. 

In the same figure we also show the result for the annealed approximation for the line domain based on Eqs. 
(|33fl - l|36ll . We see that this provides a lower energy than the quenched result for positive j3 and higher for negative /3. 

The fluctuation in the energy ((e — (e)) 2 ) at high temperature can be computed using the order (3 term in Eq. 
(|22|l . This takes different values for annealed and quenched cases, but agrees in each case with the respective values 
computed using the combinatorial approach of reference |l6j . 



In two dimensions the domains we consider are the sphere, torus and the unit square or box. 

It is convenient to treat a sphere of unit radius so the basic equations need a little modification to deal with a 
domain V whose size is not one but V. In effect, all integrals appearing in the formalism are normalized by the 
volume, and the final expression for the energy E((3,V) scales as: 




(76) 



1. 



One Dimension 



2. Two Dimensions 



(77) 



The annealed/quenched equations have a constant solution with: 




(78) 



leading to energy (normalized for a unit size domain) 

2/3 



+ 



(79) 



6a (3 2 + 4tt 
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One Dimension 



box: annealed 
box: quenched 
box: Monte Carlo 
ring: quenched=annealed 
ring: Monte Carlo 




Beta 

FIG. 4: Theoretical prediction for the average energy e for one dimensional TSP on ring and line domains as a function of 
j3 (solid line) compared with the Monte Carlo simulations (solid circles). The results of the annealed approximation for the 
line domain is also shown. Negative /3 corresponds to the Maximum TSP problem. Error bars based on 20 realizations of the 
quenched points are smaller than the symbol sizes. 



Note that in the limit (3 — » — oo, the energy becomes the half circumference corresponding, as is the case for all the 
closed domains, to the maximum distance two points can be apart. 

For a torus we return to unit size domain normalization and find that the annealed/quenched equations yield 



P 2 Jo 



tt/4 



(1 



-/3/2cos0 



f3e- 



-/3/2cos0 



2cos6» 



-)d0 



2ir 8 /• ln(1+yl) 1 + fcoshw 
WW Jo coshw f 



- cosh v \ 



dv 



In the high temperature limit the integral can be evaluated analytically leading to: 

1 



e a (/3 = 0) = - (V2 + ln(l + \/2)) = 0.382597858 



(80) 



(81) 



At very low temperature we can do a saddle point near v = in the second version of the integral. This gives 
A — ► 2tt / (3 2 and e — > 2//3. The same result comes more simply from realizing that only short distances contribute in 
the original 2 dimensional integral. A similar argument can be used for large negative [3 to give and 
e-»l/V§-2/|/3|. 

For the traditional TSP on a unit square domain, the quenched result is different from the annealed approximation 
and there is little hope of a general analytic solution. Particular values, such as the high temperature j3 = value 
may be evaluated (given patience with a 4-D integral). In a later section we derive expansions for large \(3\. 

Figure (J5J shows the average energy for each of the three domains. For the sphere this is given by l|79l) . for the torus 
it is based on numerical integration of i|80fl and for the box we resort to an iterative solution of the original quenched 
equations. The accuracy of this iterative technique is confirmed by reproducing the results for the other domains. In 
all these cases we have also performed Monte Carlo simulations and obtain excellent agreement with the theory. At 
large positive (3 the topology starts to become unimportant and each domain has energy ~ 2/(3 as expected for a two 
dimensional TSP. 
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Two Dimensions 




FIG. 5: Theoretical prediction for the average energy e for two dimensional TSP on a sphere, box and torus as a function of /3. 
Negative j3 corresponds to the Maximum TSP problem. 



D. Descent model 



In the case of the descent model, V(x) is clearly not an even function of x, and the general approach presented 
above does not apply, however the value of \ q as a functional of \x can be explicitly computed. Despite the fact that 
we know the annealed approximation produces the correct result we shall pursue our method in this case as it is 
rather instructive to do so. 

Defining R{x) = J* exp(— /i(y) /2) f^\y)dy in Eq. i|6TJ|) we find 



dR 



exp(-M)^ = V [(1 - exp(-f3))R(x) + exp(-/3)i?(l)] 



Now we define y(x) — J? dy exp(— fi(y)) to obtain 



This can be solved giving 



^ = A" 1 [(1 - eM-0))R(v) + exp(-(3)R(y(m 



!) - 1 )Srv7 + l = exp(A- 1 y( a; )(l-e X p(-/3)) 



Now setting x = 1 in the above gives 



which yields 



We thus find that 



exp(/3) = exp (A^i/(1)(1 ~ ex P (-/3))) , 



1 - exp(-/3) 
5\n{\ q ) 



dyexp(-^(y)). 



exp(-fi(x)) 



W fidye X p(-Li(y)) 



(82) 



(83) 



(84) 



(85) 



(86) 



(87) 
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FIG. 6: Predicted average energy for potential V(x) = x 2 as a function of (3 (solid line) against value measured from Monte 
Carlo simulations (filled circles) 

A solution (there is again a family related by a constant factor giving the same action) to the saddle point equation 
is 



exp(-/x(z)/2) = yjp q (x). 
This now yields 

0Fn 

= — ln(l — exp(— 0)) + ln(/3) + terms independent of (3, (89) 
This is the same result as the annealed calculation of the precedent section as expected. 

E. Other one dimensional potentials 

As an additional numerical verification of our method we have considered the potentials V(x) = x 2 and also 
V(x) = — ln(| x| ) on the line domain T> = [0, 1]. The latter potential was only considered at positive temperature as 
it is ill defined at negative temperature. The comparison of the predictions of our method against results obtained 
from Monte Carlo simulations (carried out with the same protocols as for the TSP case) are shown in Figs. © and 
@. The agreement is again excellent. 

IV. GENERAL ZERO TEMPERATURE BEHAVIOR IN ONE DIMENSION 

In some one dimensional cases we may analyze the low temperature behavior of Eq. I|66(l and extract the low 
temperature energy of the system analytically. We write s(x) = exp(— (3w(x))t(x) and thus Eq. (|66l) becomes 

exp hm*)) t(x) = J dy ^(-m^y) + My)) _ (90) 

We now assume that lri(£(x))//3 — > as /3 — > oo which permits us to evaluate the integral on the right hand side of 
Eq (|rjU|l via the saddle point method: 



exp(— (3w{x))t{x) ~ exp (— [3r{x)) I d£exp 



IK 



2 



2 (V"(x-x*(x))-w"(x*(x))) 



(91) 
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FIG. 7: Predicted average energy for potential V(x) 
Monte Carlo simulations (filled circles) 



= — ln(|ai|) as a function of /3 (solid line) against value measured from 



where 

r(x) = V(x — x*(x)) — w(x*(x)) = min {V(x — y) — w(y)\ , (92) 

1/6(0, t) 

and the point x* (x) is simply the point about which the action in the saddle point is minimal, the fluctuations being 
integrated about this point. It is at this point tempting to suggest a tentative physical interpretation of x*{x) as the 
optimal point that the polymer jumps to if its current position is x, that is to say if the monomer i is at x then the 
optimal position for monomer i + 1 is at x* (x) . 

The function w(x) is thus determined by r(x) = w(x) i.e 

w(x) — min {V(x — y) — w(y)} (93) 

2/6(04) 

Intriguing we will see that by effectively guessing some (local) heuristics for x*(x) we will be able to obtain some 
solutions to Eq. I|93(l . Putting all this together we obtain the equation for t 

^)^)) = { w j7 ¥^^r^wm ) 2 ' (94) 

this only being valid if V"(x — x*(x)) — w"(x*(x)) > on all but a set of measure zero and when, again on all but 
a set of measure zero, the minimizing point occurs within the domain [0,1]. To simplify our analysis we shift the 
domain [0, 1] to the domain [— i, — i], by symmetry we now expect that w(x) — w(—x) and t(x) = t(—x). The shifted 
equation for w is simply 

w(x) = min {V(x — y) — w(y)} (95) 



2 ' 2 > 



The energy is now given by 

d n- 

dx w(x) - 2 — / dx Mt(x)) (96) 

If indeed x*(x) is the point which is the optimal to jump to from x we expect a one to one correspondence between 
x and x* (x) in order to generate a uniform annealed distribution. If this is indeed the case, then for any function F 
on [0, 1] we will have 



(97) 
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Using Eq. $Q and Eqs. JHSJ and §52) we find 

da; wfx) — 7-— / (ix ln(t(a;)t(x*(a;))) 

1 

= 2j* i dxv>(x) + ± (98) 

Thus, given that the conditions stated above all hold, the correction to the zero temperature energy at low tempera- 
tures takes a remarkably universal form. 

We first consider the case where V is a purely attractive potential with a minimum at x = 0. Taking the idea that 
x*(x) is the optimal jump from the point x we expect x*(x) = x. This will imply from Eq. (|93|) that w(x) — V(0)/2. 
This solution can be seen to work when plugged back into £'a. (l93f) when V'(0) = and V"(0) > 0. We thus find from 
Eq. fljjBJ that 

e ~ v{0) + ip (99) 



We now consider the case where the potential is everywhere repulsive. On the interval [— ^, ^] the greedy algorithm 



I - 1 1 

described earlier amounts to making the choice x*{x) = —x. This choice implies that 

10(35) = V(2x) - w(x) (100) 

and for the solution to be valid we must have that 

-V' {2x) - w'(-x) = 0. (101) 

If w(x) — w(—x) then we have w'(x) — —w'(—x) and so the above implies that w{x) = V(2x)/2. For this solution to 
be valid we must have that V"(2x) < and thus it only holds for concave potentials. 
The low temperature energy is thus given by 

e = JdxV(x) + ± (102) 
The ground state energy is clearly that given by the greedy algorithm 

e GA = I das V(x). (103) 



Another heuristic for finding the maximal path in the case of repulsive potentials is to take a jump size of constant 
size A. When there are no minima the best value for A is 1/2, this strategy clearly the minimizes the energy at each 
jump subject to the constraint that it must be possible from any position x. One simply adds a very small noise r\ 
to each jump of 1/2 to generate the required uniform distribution of monomers on [0, 1]. We call this the half jump 
algorithm and it is clearly not at all greedy ! We thus take x*(x) so that \x — x*(x)\ = 1/2. When x > this implies 
x* (x) = x — 1/2 and hence 

w{x) = V{ 1 -)-w{x- 1 -) (104) 
The explicit solution to this equation is w(x) = a\x\ + b, where 



b=V{\)-~b, (105) 



and the condition to have a minimum implies that 



V'(^)-a = 0, (106) 



and V" (1/2) > 0. Thus the function V cannot be concave near x — 1/2 and 



(107) 
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FIG. 8: Value of — \n(s(x))//3 obtained by numerically solving Eq. I|660 at /3 = 30 (solid line) for the potential V(x) 
Also is shown the zero temperature prediction for this function (dashed line). 



As expected this solution gives 



(108) 



which is obviously the ground state energy given by the half-jump algorithm energy. 

A potential where the above solution is possible is V(x) = — ln(|x|). Numerical solution of Eq. I|tj6f> & t l° w 
temperatures converges to the solution found above. The predicted value of the ground state energy is eqs — m (2), 
this value is compatible with the Monte Carlo simulations for this potential shown in Fig. ||7J|. Clearly the greedy 
algorithm is a bad strategy for the potential V(x) — — ln(|x|), this is because at the end it must link points very close 
to each other situated near x = 1/2, thus giving a very large contribution to the energy at the end of the algorithm. 
The greedy algorithm gives an energy e = — Jdxln(x) = 1, which is indeed higher than the ground state we predict 
analytically and not compatible with our Monte Carlo simulations. When V"(|x|) > everywhere in [0, 1], Jensen's 
inequality tells us that (V(X)) > V(X) for X distributed on [0, 1]; when this distribution is uniform this implies that 
tQA > £ha and hence the half jump algorithm is the most efficient. In the case where the potential is concave the 
greedy algorithm is the most efficient. 

We note that the case of the Maximum TSP is an intermediate case where V"(x) = and in this case to a — £ha 
and the forms of u(x) in these two cases coincide. As a check of this asymptotic analysis we have numerically solved 
Eq. (jnnjl for the potential V(x) = -\x\ at (3 = 30. Shown in Fig. © is w*(x) = -ln(s(x))//3 where s(x) is the 
numerical solution of Eq. i|66|) at (3 = 30. Also shown is the predicted zero temperature limit of w*. The agreement 
improves on increasing /3 but limitations of numerical accuracy are attained if (3 is taken to be too large. 

The finite temperature corrections for both the TSP and Maximum TSP are different to those of the cases considered 
thus far as V"(x) = for x ^ 0. At positive temperature it is clear that w(x) = 0, as the saddle point is at x*(x) = x. 
Here, because V"(x) — 0, we do not expand the the terms in the exponential of the integral about y = but we do 
carry out the expansion of the term l/t(y), we thus write 



t 2 (x) w / dx exp (—/91a; — y\) = \ ( 2 — exp 



'(-f - P x ) ~ ex P(~^ + fa] 



(109) 



As (3 — > oo we have 



dx ln(s(x)) 



o 
1 

w 



dx In ( 1 — — exp(— j3x) 



HI 



HI 



(110) 
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This yields 



y In (2) . (Ill) 



for the TSP as f3 — > oo. We have checked that this result agrees with the asymptotics of our numerical solutions. 

The analysis for the Maximum TSP is more involved. Differentiating Eq. 165|) twice we see that the function s(x) 
obeys 



s"(x) = /3 2 s(x) + -^. (112) 
s{x) 

We now make the substitution s(x) = exp(/3\x\)t(x) to find 

t"(x) + 20 S gn(x)t'(x) = -2j35(x)t(0) + 2 P e M~W\A) _ (n3) 

t{x) 

Now in the limit /J^oowe have 2/3exp(— 2/3|x|) « 8(x). For large (3 we thus have 

*"(*) + 2fcgn{x)t'{x) = S(x) (-2(3t(0) + ^-J . (114) 

The above has solution t{x) = A + B exp(— 2/3|ir|) and the jump conditions at the origin yield the relation 

A 2 -B 2 = ± (115) 

An extra relation between A and B is found from examining the integral equation for t at x = 0, which in this limit 
gives: 

A + B = ± (116) 
and in the limit of large (3 we find A ^ B ^ 1/ \[2 and we find that the large [3 behavior of the energy is 

1 7T 2 , „ s 

^-2 + W 

To end our analysis of the low temperature limit we will consider the Maximum TSP in higher (d) dimensions, 
specifically on the hypercube [0, l] d . Consider the generalization of the greedy heuristic. Here we start on the 
outermost layer of points in the hypercube and we join points on this surface to those that are diametrically opposed. 
The procedure is then repeated eating away the hypercube until we arrive at the center. Shifting the domain to 
[— as before, this entails matching the point x with —x. This generalized heuristic was shown to give the 
optimal path length for d = 2 [23] • The ground state energy generated by this generalized greedy heuristic is clearly 



e GA = 2 d [ dx V(2x) 



= ~2 d+1 f dx |4 (118) 

This general formula gives cqa = —0.5, , —0.765196, , —0.960592 in one, two and three dimensions respectively. The 
solution w(x) — —x is in fact a solution to Eq. l(9"3|l and thus gives these ground state energies. This can be easily 
verified as we note that the function 

h{x) = \y\-\x-y\ (119) 

is bounded as 

h(x) > -\x\ (120) 

by the triangle inequality. The bound is achieved at y — —x, confirming that w(x) = —x is indeed a solution. We 
note that this solution exists in any domain T> (centered at the origin) satisfying the property that if x € T> then 

-x e v. 



V. CONCLUSIONS 



We have discussed the statistical mechanics of models whose phase space is the set of permutations of N objects 
characterized by quenched positions r.; . The Hamiltonians are functions of the neighboring elements in the sequence, 
and thus a given sequence can be interpreted as the energy of a polymer ring or closed random walk which visits 
all points in the quenched distribution once. We analyzed the cases corresponding to several well studied problems 
including the traveling salesman problem the descent problem and the quenched Rouse model. 

The annealed approximation was first considered and illustrated for some one-dimensional cases. For the TSP on a 
ring and the descent problem this annealed approximation gives the correct quenched result. For the descent model, 
this is because the effective potential between neighboring monomers, when the system is viewed as a polymer with 
interactions between consecutive monomers, is scale free and independent of the quenched distribution of the random 
points. For the TSP on a ring, the reason is less clear, but agrees with expectations from replica studies of the 
independent link approximation and continues to hold for symmetric closed domains in higher dimensions. However 
in general we expect the annealed approximation to fail at all but infinite temperature. This is because the points 
are allowed to evolve dynamically to lower the free energy of the system and the resulting thermodynamic distribution 
will not be the same as their original quenched distribution. 

We then showed how the quenched calculation could be carried out and confirmed its predictions for both one and 
two dimensional TSP examples with Monte Carlo simulations. Physically the method we introduced corresponds to 
imposing a fictitious site dependent chemical potential on the distribution of a set of dynamical variables in the 
presence of the original interaction Hamiltonian. This chemical potential is then chosen to ensure that the annealed 
distribution of the positions of these dynamical r^, denoted in this paper by p a (r), is the same as the quenched 

distribution of the quenched random variables denoted by Pq(i'). The method is exact in the thermodynamic 
limit (corresponding to high density where the length of the interval is held constant) for any quenched distribution 
p q (r) and interaction potential V(t). O n of the most intriguing observations made here is that in the annealed 
approximation wc arc lead to consider a linear eigenvalue problem to solve the thermodynamics as we use a transfer 
matrix approach, however in the quenched calculation we are lead to consider a non-linear integral equation. Although 
we managed to avoid the sometimes rather opaque replica method in our treatment, it would be interesting to see if 
the appearance of the non-linear eigenvalue equation could be interpreted or re-derived within the replica formalism. 
The results of our calculations were then confirmed by comparing them with Monte Carlo simulations in a variety 
of models. In the cases we have considered so far we have seen no evidence for any phase transition on lowering the 
system's temperature. It is possible however that in higher dimensions and with certain interaction potentials V that 
a phase transition does occur. We recall that the a directed polymer in dimensions greater than two exhibits a finite 
temperature phase transition pi) . 

Particular attention was paid to the zero temperature limit where Eq. I|93|l needs to be solved. A number of 
solutions were found which, although we did not prove uniqueness, are compatible with our numerical simulations 
and also a rigorous result for the Maximum TSP in two dimensions. The solution of Eq. I|93|) was carried out by 
trying out different heuristics to construct the optimal path, it is possible that in more complex situations the method 
used here could serve as a useful indicator for such constructions. 

Finally the idea of treating quenched variables as effectively annealed variables and then adjusting their Boltzmann 
weight in order to recover self consistently the original quenched distribution may prove useful, either as an exact 
or approximate method in other problems involving quenched disorder. Indeed Morita's pioneering work used this 
idea in an approximate context, here we have shown that the procedure can be carried out exactly for this type of 
permutation based combinatorial optimization problem, as is also the case for some one dimensional spin models |l9j . 
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